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We report on a detailed examination of numerical results and analytical calculations devoted to 
a study of two holes doped into a two-dimensional, square lattice described by the t — J model. Our 
exact diagonalization numerical results represent the first solution of the exact ground state of 2 
holes in a 32-site lattice. Using this wave function, we have calculated several important correlation 
functions, notably the electron momentum distribution function and the hole-hole spatial correlation 
function. Further, by studying similar quantities on smaller lattices, we have managed to perform a 
finite-size scaling analysis. We have augmented this work by endeavouring to compare these results 
to the predictions of analytical work for two holes moving in an infinite lattice. This analysis relies on 
the canonical transformation approach formulated recently for the t—J model. From this comparison 
we find excellent correspondence between our numerical data and our analytical calculations. We 
believe that this agreement is an important step helping to justify the quasiparticle Hamiltonian, and 
in particular, the quasiparticle interactions, that result from the canonical transformation approach. 
Also, the analytical work allows us to critique the finite-size scaling ansatzes used in our analysis of 
the numerical data. One important feature that we can infer from this successful comparison involves 
the role of higher harmonics in the two-particle, d-wave symmetry bound state — the conventional 
(cos(k x ) — cos(ky)) term is only one of many important contributions to the d-wave symmetry pair 
wave function. 



I. INTRODUCTION 

The behaviour of mobile holes in an antiferromagnetic (AF) spin background has been a subject of intensive study, in 
part because of its possible connection to high temperature superconductivity. The ubiquitous structural components 
of such materials are the Cu02 planes, and the description of carriers in these planes is the theoretical focus of this 
paper. We consider the so-called t — J model MM, for which the holes correspond to the Zhang-Rice singlets [pj, 
mobile vacancies created by doping a single Cu02 plane. A microscopic representation of this model is 

Ht-j = -tj2 (4^ + H - c -) + J E( S * • s ^ - i n ^)' f 1 ) 

where (ij) denotes nearest neighbour sites, and c| CT , Ci a are the constrained operators, cur = Ci a (l — c\ _ .c, ) _ CT ). The 
■ ratio of the AF exchange constant to the hopping amplitude is believed to be J/t ~ 0.3. 

Aided by recent angle-resolved photoemission experiments followed by extensive comparisons between theory 
and experiment J^] , it is now recognized that this simple, near-neighbour hopping Hamiltonian on its own is insufficient 
to fully represent the true single hole state of the real Cu02 plane. Hoppings between more distant neighbours are 
required as are more complicated three-site spin-dependent hoppings. Possibly, the full three-band microscopic 

Hamiltonian is necessary 

Despite the potential inadequacy of this Hamiltonian to represent completely the Cu02 planes, it is still the simplest 
model that captures the important antiferromagnetic correlations of a weakly doped antiferromagnet Thus, it is 
crucial that the properties of this model when doped are elucidated. 

The Hamiltonian in Eq. (|l|) has been investigated intensively by different analytical and numerical methods, and 
we believe that a consistent picture has emerged from these studies. Some results of the exact diagonalization on 
small clusters, such as the energy spectrum and quasiparticle residue, have been found in the excellent agreement 
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with the results predicted by advanced analytical theory JlO(]. In contrast with this, a large amount of the numerical 
data for the correlation functions are not well understood or require further explanation. Such an explanation, when 
completed, could help to define (or justify) the correct quasiparticle model for the system of strongly interacting holes 
and spins at low energies. It is in this manner that we unite our analytical and numerical work in this paper. 

A common and apparently successful theory of a single hole moving in an AF aligned background involves the 
so-called spin polaron Jll| jlj] . According to the spin-polaron idea, the hole in its movement disturbs the magnetic 
background that one can formally describe as the strong coupling of the hole and spin degrees of freedom. This makes 
this problem similar to the well known strong coupling electron-phonon one. However, in spite of the qualitative 
similarity of these two polarons, there is an essential difference between them. If the phonon polaron can be considered 
as an almost static object of the shifted ions with the electron in the centre, the spin polaron is the "spin-bag" with the 
moving hole inside. One of the statements of the present paper is that this feature of the spin polaron is responsible for 
the absence of the direct similarity between the answers which theory provides for quasiparticles and the numerically 
obtained data for bare holes. A similar conclusion, using a different analytical approach and numerical results for 
smaller clusters, was reached by Eder et al [fL5| P~T|| - Later, similar remarks were made by Riera and Dagotto [Q. 

In this paper we combine analytical and exact diagonalization (ED) numerical results of the one- and two-hole 
problem to provide a comprehensive study of these important systems. Computationally we have managed, for the 
first time, to determine the two-hole ground state for two holes doped into the 32-site, square cluster used by two of 
us in a previously published numerical work JTof . We find that the lowest energy state is a P = bound state with 
d x 2_ y 2 symmetry. We have characterized the ground state by evaluating a number of important expectation values, 
notably the electron momentum distribution function (EMDF), and the spatial pair correlation function. 

We have found that an effective quasiparticle Hamiltonian, originally proposed by Belinicher, one of us, and Shubin 
p9| , may be used to calculate the same expectation values that were obtained numerically via ED. Further, these 
quantities are remarkably similar to those obtained via ED. This gives strong support to the appropriateness of this 
quasiparticle Hamiltonian. 

Previous analytical work on the low-energy physics of the two-hole system generally describes it in terms of mod- 
erately interacting spin polarons. This analytical work shows that the dominant effective interactions between spin 
polarons come from the short-range nearest-neighbour static attraction and spin-wave exchange, the latter leading to 
a long-range dipolar interaction. These interactions are attractive for d-wave states and strongly repulsive for s-wave 
states. The purpose of this paper is to use the ED results to provide support for this description of the internal 
structure of the quasiparticles, and indirectly for the above-mentioned description |l9f| of their interactions. 

We will first describe the present status of the t — J model studies in §11. Section III discusses in detail the numerical 
data available for the ground state correlation functions. Then, §IV summarizes briefly the analytical results of 
relevant previous work and displays the details of the present calculations. Section V focuses on the comparison of 
the analytical and numerical results, and in §VI we present our conclusions. 



II. PREVIOUS t - J MODEL STUDIES: 



A. Analytical results 

Theoretical studies of the t—J model have resulted in a clear understanding of the nature of the low-energy 
excitations for the system near half filling. The charge carrier created by a hole introduced in an AF background 
is described as a spin polaron, viz. as a quasiparticle consisting of a hole and a cloud of spin excitations. The AF 
spin-polaron concept was put forward in earlier works of Bulaevskii et al [ pT| and Brinkman et al p(i| ], and then 
developed in a number of more recent papers ||l^ [l^ , ^l| "^8| using different techniques. 

The main conclusion of these papers was that the spin polaron in an AF background is a well defined quasiparticle 
with a nonzero residue and a specific dispersion law. The dressing of the hole leads to the narrow quasiparticle 
band with a bandwidth ~ 2 J for realistic J/t, band minima at k = ±(7r/2, ±7r/2), and a heavy effective mass along 
magnetic Brillouin zone (MBZ) boundary. Most of these features of the spin polaron were found to be robust under 
generalizations of the t — J model p9[ , including further neighbour and three-site hoppings, and for a wide range of 
J/t ratio. 

The single- hole problem has been treated analytically in detail [pl]- pTpO|j22] - p7| , pO| . Grouping these efforts, two 
approaches in treating this problem were used: (i) the self-consistent Born approximation (SCBA) (e.g., see Refs. 
|L3| . [r4|j22| ) , and (ii) the so-called "string" approach (e.g., see Refs. |Il| , ^o| , p^ , p5| -p7| ) . Recently, a relationship between 
these two has been established [^lj. We briefly review these methods with an eye to understanding how well they 
might be able to describe the two or multi-hole problem. 
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The SCBA method utilizes a property of the hole-magnon interaction, namely the absence of the lowest order 
correction to the Born approximation series for the single-particle Green's function Jl3],[l4|,[22| . The attractive feature 
of the SCBA approach is that the essentially exact single-hole spectral functions can be evaluated quite easily using 
simple numerical calculations. Recently, the detailed structure of the single- hole ground state and behaviour of 
different correlators has been studied using SCBA p2] . Unfortunately, already in two- or many-hole problems much 
more involved numerical and analytical efforts are required |33f| . 

The string approach is based on the idea that in an AF background a hole will be confined by an effective potential 
created by overturned spins ("strings"). Formally, the real-space variational ansatz for the polaron's wave function 
with the strings of different length is considered to reproduce this tendency. In spite of the considerable success of the 
string approach for the single- hole problem |25|] , there are some problems which make the use of it as a candidate for 
a quasiparticle theory for the t — J model questionable. First of all, the string method uses the real-space approach 
which does not treat properly the long-range dynamics of the system. The next problem is that the method starts 
from the Ising background and includes fluctuations on a perturbative basis, whereas the fluctuations are strong in 2D 
and must be "build in" to the real ground state of the spin system. The third problem concerns the necessity of the 
normal anticommutation relations of the quasiparticle operators. If hole is a fermion a unitary transformation, which 
diagonalizes the Hamiltonian and dresses the hole by the spin excitations, would not change its statistics. Within the 
variational "string" and other approaches, one works with the wave functions and usually identifies the wave function 
of the quasiparticle with the operator of the quasiparticle. This leads to the absence of the commutation relations for 
these operators and to troubles with the proper normalization and orthogonality of the states already for the two-hole 
problem |l8|,Q. Because of that one cannot correctly derive the effective polaron-polaron interaction term using a 
single-hole wave function. 

These difficulties notwithstanding, several attempts to address the two-hole ground state have been made. Work 
based on the string approach have led to some qualitative understanding of the problem |l7],^3|,[35| . 

The investigation of the interactions between quasiparticles in the t — J model is a subject of prime interest in 
the context of magnetic pairing mechanism. However, studies of this problem show much less convergence than the 
single-hole problem. Other work involves the formulation of an effective model for spin polarons propagating in the 
AF background |3(]|36| p9fl . From the RPA treatment of the Hubbard model in the strong-coupling limit, the model 
of "spin-bags" interacting via longitudinal magnetization fluctuations has been proposed J3lj. A phenomenological 
model for the vacancies coupled by the long-range dipolar twist of the spin background has been also worked out 
po| , [fOf using a semi-classical hydrodynamic approach. Inspired by the numerical evidence of the hole-hole d-wave 
bound state and the Van Hove singularity in the single-hole spectrum, taking the simplest phenomenological form 
of the interaction, an AF Van Hove model has been put forward ]39| . Using an ansatz for the quasiparticle wave 
function [^7| , the "contact" hole-hole and the residual hole-magnon interactions have been obtained |4l]j42[ | , and then 
the effective Hamiltonian for the polarons and long-range spin- waves has been presented [p6| . 

Even though most of these theories were formulated on a phenomenological and semi-phenomenological basis, they 
provided two key interactions leading to pairing in the t — J model. One of them is the effective hole-hole static 
attraction due to minimization of the number of broken bonds found from placing two holes at nearest-neighbour 
sites (sometimes referred to as the "sharing common link effect"). The other is due to spin- wave exchange, and leads 
to a dipolar-type interaction between holes (3^j3^,Q . 

Quite recently, a new approach to the derivation of a quasiparticle model from the t-J model has been developed 
ipTsfl . It used a generalization of the canonical transformation (CT) approach of the Lang-Firsov type. An effective 
Hamiltonian for the spin polarons includes in itself both types of the hole-hole interactions mentioned above in 
a natural way. Some details of this approach are presented in Section IV. In Ref. results for the single-hole 
properties have been compared with ones of the SCBA calculations and an impressive agreement has been found. 
This is supported further by the idea that the "canonically transformed" quasiparticles are close to exact t — J model 
ones. Even though CT approach is less controlled than the SCBA one, it solves naturally all the problems mentioned 
above and allows one to derive the quasiparticle Hamiltonian for interacting spin-polarons from the original t-J model. 
Thus, in this paper we compare the predictions obtained from this Hamiltonian to our numerics, and in this way we 
critique the description of the interactions between quasiparticles that follows from the CT approach. 



B. Numerical studies 



ED studies of the t — J model doped away from half-filling on small clusters with periodic boundary conditions 
are an important source of unbiased information on the low-energy physics of this system. One- and two-hole ground 
states have been investigated in great detail on the 16- (4 x 4) 18- (VT8 x \/l8), 20- (V20 x a/20) [g3^2j, 

and 26-site (a/26 x \/26) ]63|-|6l| clusters. Although some of these results converge, at least partially, these clusters 
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suffer from various finite size problems. The 20- and 26-site clusters do not have the full rotational symmetry of the 
square lattice. Therefore, they do not possess important reciprocal lattice points along the high symmetry directions 
in the first Brillouin zone. For example, the important reciprocal lattice point (ir/2,w/2) does not exist in the first 
Brillouin zone of the 18-, 20-, and 26-site clusters. This causes the ground state momenta of the one-hole state to 
be different from the predicted ±(7r/2, ±7t/2) points. Although the 16-site cluster has the (tt/2, 7r/2) point, it has an 
additional symmetry which causes an accidental degeneracy of the levels at (ir/2,ir/2) and (it, 0) for one hole, and 
between (0,0) and (tt, 0) for the two-hole problem |52]. Attempts to remedy the missing (7r/2,7r/2) point have been 
made by using the non-square 16-site (a/8 x a/32) Jnf and 24-site (a/18 x a/32) @ clusters. 

Previous results on the single-hole problem show that the quasiparticle peak at the bottom of the spectral function 
can be expected to survive in the thermodynamic limit |IoU4lf , fi^ , |64|| . The corresponding quasiparticle band is narrow 
(of the order of 2 J in the "physical" region J/t < 1) and the band minima are shifted to the ABZ boundary. However, 
due to the previously mentioned deficiencies of the 16-, 18-, 20-, and 26-site clusters, none of them can actually show 
that the quasiparticle minima are at the ±(tt/2, ±7r/2) points; these wave vectors are those predicted by numerous 
theoretical studies 

The smallest cluster which has the full rotational symmetry of the square lattice, contains the (tt/2,tt/2) point, 
and is free from the spurious degeneracies mentioned above, is the 32-site cluster (a/32 x a/32) — see Fig.|l| Also, it 
is the largest such system which can be solved using modern computers. These calculations involve finding the lowest 
energy states of matrices with dimensions of up to 300 million. 

Recently, some results for the single-hole problem have been published for this cluster by two of us fTo) . These 
numerics showed that the effective mass around the minima is anisotropic, and that the quasiparticle residue is rather 
small for realistic J/t, both in excellent agreement with analytical predictions. Further, the full dispersion relation 
predicted by analytical work based on the SCBA jll].|l4|,|22j is found to be in excellent agreement with ED numerics 
on this cluster jlQl. 




FIG. 1. The 32-site cluster in (a) direct and (b) reciprocal space. Empty circles are included to mimic the periodic boundary 
conditions used in our studies. The solid lines in (b) show the important directions in k-space displaying the high symmetry 
of the cluster. The dashed line in (b) borders the first magnetic Brillouin zone. 



ED results for two holes on finite clusters consistently show that they are coupled in a bound state with d x 2_ y 2 
symmetry in a wide range of J/t @,|^,|^,|§j50|j5|,|^,|^,|l]J^|9|,^2|, in agreement with earlier ED data for the 
16-site Hubbard model J73| , [7^ |, the Green's Function Monte Carlo studies on 8 x 8 cluster JT5[ , and some theories 
JTrl . pi) 36 3^j76|. Low lying states with other symmetries, as well as k ^ (0,0) states, have also been studied p9| , |6l| . 
Attempts have been made to extrapolate the binding energy to the thermodynamic limit and thus to estimate the 
critical value of J/t for the formation of a bound state |6^^]. Also, some knowledge concerning the internal structure 
of the bound state is known Bq,|7q]. Lastly, the electron momentum distribution function has been investigated in 
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order to search for Fermi-like discontinuities. Drastic differences between the single-hole and two-hole cases have been 
noted 51 5^J|^]. Finite size and J/t scalings of this quantity have also been studied 58 59|j6^]. 

Nevertheless, the above-mentioned results are of limited usefulness simply because of the systematic error introduced 
by the lower symmetry of the clusters with 16, 18, 20, and 26 sites. Clearly, the 32-site cluster would augment such 
studies. Further, with the collection of all such clusters, and some analytical guidance regarding the correct scaling 
laws, information on the thermodynamic limit, viz. a density of zero (two holes in an infinite 2D square lattice), 
would be accessible. 

Section III will summarize the results of the ED studies of the single- and two-hole problems that we have obtained 
on the 32-site cluster. Some of the single-hole results have been published previously |l^], and we only mention those 
results that are crucial to our scaling analyses. A brief summary of a portion of the two- hole comparison to the CT 
Hamiltonian was presented in Ref. p% . 



III. NUMERICAL WORK 



Most of the results presented in this section are obtained by ED on the 32-site cluster with periodic boundary 
conditions for the realistic value J/t — 0.3. Results at different J/t, as well as on smaller clusters published previously 
|3l],[37|], are also used in the discussion of the finite size scaling (FSS), bound state energies, and correlation functions. 

To look for evidence of hole binding in the low energy states, we calculate the two-hole binding energy = 
E2 — 2Ei + Eo, where E\ and Eq are the ground state energies with one and no hole respectively, and E2 is the 
energy of the two-hole state. Further, two expectation values that we are interested in are defined as follows: (i) 
The electron momentum distribution function (EMDF) is given by (u^) = (cj^Cko-), where , Ck CT are the Fourier 
transform of the constrained operators, (ii) The spatial distribution of holes in the bound state is characterized by 
the pair correlation function defined as 

C{r) = N N (r) " ni){1 "*M*-'V>> ( 2 ) 

where Nh is the number of holes, and Ne(t) is the number of equivalent sites at a distance r from any given site. 

Before presenting the FSS of the EMDF, we discuss what kind of finite size behaviour one can expect. The EMDF 
is expected to show how hole doping changes the uniform value of (riko-) = \ obtained in the half-filled case. In a 
system of free particles, a hole with a certain momentum will manifest itself as the complete suppression of (fik<r) to 
zero at this momentum. In systems with interaction whose physics can be described in terms of the quasiparticles, 
this suppression will be proportional to the quasiparticle residue Z^, and is almost independent of the cluster size; 
the rest of the hole weight will be distributed among the other available k-points. Consequently, the more k-points a 
system possesses the less hole weight each k-point will carry. Therefore, in the single-hole problem we anticipate that 
(n-ka) for the ground state momentum P to be suppressed by a constant proportional to Zp. Further, we expect that 
the deviation from the half-filled value, 

(Sn klJ ) = (n ka ) - i, (3) 

will scale as 1/N at all other points until (roughly) the physics of the system does not change with the cluster size, 
i.e., when size of the quasiparticle is smaller than the cluster size. 

The same argument should apply to the bound states of the two-hole problem, where we predict (Sn^) to scale as 
1/N at all k-points. 



A. Single-hole case 



We wish to provide a FSS analysis of certain quantities for the two-hole ground state. To this end, we present new 
results for the one-hole problem that will facilitate such work. 

Figure shows the EMDF of the single-hole ground state on the 32-site cluster at J/t = 0.3, which has total spin 
Sf ot = +5 and momentum P = (7r/2,7r/2). Due to the non-zero momentum of this state, the only symmetry its 
EMDF has is a reflection about the "main diagonal" ((— 7T, — it) «-» (n, it) line). 

Several features of the EMDF are worth noticing. First, (rikx) nas a "dip" at the GS momentum P. Earlier studies 
of the J/t dependence of the intensity of this "dip" have left no doubt about its direct relation to the quasiparticle 
weight Zp p8| . Second, the EMDF deviates significantly from its half- filled value for both spin directions across the 
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FIG. 2. The EMDF for the single-hole ground state having momentum P = (w/2,ir/2) and 5f ot = +1/2 at J/t = 0.3. 
The numbers are the electron filling factors (a) (n^f), (b) (rikj.} at the corresponding k-points. "Antidips" at (— 7r/2, — 7r/2) 
described in the text are denoted by the dashed boxes, and the "dip" in (b) at the ground-state momentum is highlighted by 
the solid box. 

entire Brillouin zone. This background has a maximum at (0,0) and a minimum at (7r,7r), which is a biproduct of 
minimizing the kinetic energy of the system Jlfj| ] . Although this "dome" shape resembles the "large Fermi surface" 
in a system of free electrons, it has different physics behind it. The discussion of this behaviour will be given in §IV. 
One observes that this dome structure in (nk|) is shifted upwards from its half- filled value, and that in (rtk|) it is 
shifted downwards ((«(o.o)t) — ( n {o.o)i) — 0.03). The difference between the maximum and minimum, 

An a = (n (0) o)a) - ("-(tt,^), (4) 

is slightly larger for a =f than for J, (An^^ ~ 0.07(0.06) at J/t = 0.3). Also An^ has a stronger J/t dependence. 
Finally, the EMDFs of both a =|, { have "antidips" at (— 7r/2, —tt/2) . They were observed earlier but no successful 
explanation has been presented. The fact that the "antidips" are always at P— Qaf supports the idea that their physics 
is somehow related to the long-range AF fluctuations in the system |]58|l . Subtraction of the "normal" background 
from (^(-Tr /2.-7t/2)o-) shows that the depth of the "antidip", 

ArianU.a = («(_. 7r/2,-7r/2)o-) _ ("(7r/2,— rr/2) CT )i (5) 

is larger for f (An a „ ti jQ) ~ 0.11(0.08) at J/t — 0.3), and An an ti,i nas stronger J/t dependence. 

In Figs. H and |we plot (fco-k) vs. the "inverse volume" of the system, 1/N, at J/t = 0.3 for all k-points (except 
P and P — Qaf) available on more than one cluster. One can see the almost perfect 1/JV scaling at all these points, 
in agreement with our expectation. Figure || shows the same plot for the EMDF at the ground state momentum. 
Extrapolation to the thermodynamic limit shows that the dip, which we expect to be Zp/2 (the factor one-half is 
from the proper normalization of the wave function), is about 0.14, or Z-p ~ 0.28. This agrees well with SCBA result, 
gSCBA _ o.271 [^3|. There is no simple scaling for the "antidips" of |(^np_Q cr )| because of the long-range physics 
involved. According to the discussion in Sec. IV they are combinations of terms with different scaling behaviours. 



G 



J/t=0.3 



J/t=0.3 





FIG. 3. 1/N scaling of the (a) |(5n kT )|, (b) |(5rakj.}| (see Eq. (§)) for k = (ir,ir) (filled circles), (0,0) (empty circles), and 
(•7T, 0) (filled diamonds) for the single-hole ground state at J/t = 0.3. The data from 16-, 20-, 26-, and 32-site clusters ((7r,7r), 
(0,0) points) and from 16-, 20-, and 32-site clusters ((n, 0) point), where these points are available, are used. 
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Cn/2.ji) 
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FIG. 4. The same as in Fig. |^ for k = (7r/2, 7r) (filled circles), (ir/2, — it) (empty circles), (7r/2,0) (filled diamonds), and 
(7r/2, — 7r/2) (empty diamonds). These points are available from 16- and 32-site clusters only. 
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FIG. 5. 1/N scaling for the |(<5rik|)| (open circles) and (C + a/N) scaling for the "dip" in [(<5rikx}| (filled circles) at the 
ground-state momentum k = P, J/t = 0.3. These points are available for the 16- and 32-site clusters only. 

B. Two- hole case 

The only zero total momentum bound state that we have found in the zero magnetization channel is a singlet and 
it has d x 2_ y 2 symmetry. Figure ^(a) shows the J/t dependence of the binding energy E\, on the 16-, 26-, and 32-site 
clusters. One can see that the absolute value of the binding energy gets smaller as the size of the cluster grows, 
and that in some region of J/t the binding energy becomes positive. Such behaviour seems to be natural in the 
presence of short-range attraction between the holes. In this case holes on larger clusters lower their kinetic energy 
due to derealization and make the bound state shallower, whereas on smaller ones they are not allowed to move 
farther apart. Further, holes on smaller clusters are forced to be in the region of the mutual attraction. Since these 
short-range interactions are believed to be of magnetic origin, the interaction energy has to scale as J. Consequently, 
the increasing importance of the kinetic energy at small J/t tends to destroy the bound state. This line of thinking 
leads to a discussion of whether or not the critical threshold of J/t for bound state formation is above or below the 
"realistic" value of J/t for the cuprates. Attempts have been made to estimate the thermodynamic limit of (J/t)\ c 
through FSS of the binding energy |67]|39|]. If we follow the same approach, we obtain the scaling shown in Fig. g(b), 
and this data shows the FSS at two representative J/t values. The thermodynamic limit of Ej, is negative at the 
larger J/t (smaller size of the bound state, larger role of the short-range interaction) and positive at the smaller J/t 
(no bound state). 

In our opinion, this approach is problematic because for at least two reasons. First, there is another important 
hole-hole interaction, viz. magnon exchange, which also leads to pairing. In fact, it is this interaction which selects the 
d-wave symmetry state. It has been shown analytically ]l9| , |30| , ^6] | that this interaction alone leads to the formation of 
a shallow long-range bound state which does not have a critical value of J/t because the interaction strength grows 
with t. Therefore, one can expects that further increase in the cluster size will not only lower the kinetic energy of 
the holes, but will also provide more sites for the holes to take advantage of the long-range attraction. The second 
reason is the absence of the evident scaling law for the binding energy. Regarding the complexity of the interactions 
involved and the tendency of the magnetic subsystem towards AF long-range order, we expect different contributions 
to the FSS of Ef, which are of different order in 1/N and of comparable magnitudes. This is especially true at smaller 
J/t when the size of the bound state is comparable to or larger than the cluster size. 
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FIG. 6. (a) The J/t dependence of the binding energy Et, in units of t, from ED studies on the 16- (diamonds), 26- (squares), 
and 32-site (circles) clusters, (b) The binding energy vs. 1/N for two representative J/t values, J/t — 0.3 (upper) and J/t — 1.0 
(lower). The solid and dashed lines are 1/N scaling using 16- and 32-site data only (filled circles) and data from all available 
clusters (open and filled circles), respectively. 



Another important quantity which shows further evidence of the hole-hole attraction in an AF background is 
the hole- hole correlation function C(r), Eq. (0). It has been studied in detail on smaller systems [^4] , p r Q|j6^j69| , |7^| . 
Figures |7|(a) and ^(b) show the 32-site ED results for C(r) at J/t — 0.3 and J/t = 0.8, respectively. In a wide 
region of J/t the strongest correlation is at the y2 distance, while the nearest-neighbour correlation is also strong. 
A density-matrix renormalization group study ]78| has also found similar physics. At larger J/t (Fig. ^|(b)) the size 
of the bound state is small: the nearest neighbour and v2 distances accumulate about 80% of the holes. However, 
at J/t = 0.3 the probabilities of finding the holes at y/B and \/2 distances are almost the same, and only 46% of the 
holes are located at the nearest neighbour and \J2 distances. The correlation decays slowly with distance at small 
J/t. Hence in the J/t = 0.3 bound state one can expect C(r) to have a longer "tail" in the thermodynamic limit. 



0.08 




FIG. 7. The spatial correlation function, C(r), for two holes doped into a square lattice described by the t-J model, for (a) 
J/t = 0.3, and (b) J/t — 1.0. Our ED results (open circles), the analytical results for an infinite square lattice (open diamonds), 
and the analytical results mapped onto a 32-site square lattice (filled diamonds), are all shown. The lines are a guides to the 
eye. In (b), analytical results for the cluster are very close to ones for the bulk, and hence are not shown. 
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FIG. 8. The EMDF for the two-hole ground state, P = 0, S? ot = at (a) J/t = 0.3, (b) J/t = 1.0 within the first quadrant 
of the BZ. 

The next correlation function which can be used to extract information on the bound state is the EMDF. Fig- 
ures ||(a),(b) show the EMDF at J/t — 0.3 and J/t — 1.0 in the first quadrant of the Brillouin zone. Since the total 
momentum of the system P is zero, the EMDF possesses the full square symmetry. Moreover, since the ground state 
is a singlet, (n^) = (riki) = (rik)- Another noticeable difference from the single-hole EMDF is the absence of "dip" 
at any k-point. This is not surprising because one would not expect the holes in the bound state to have a certain 
momentum. They will be spread over all k-points especially if the bound state is well localized in real space. 

Some of the features of the EMDF are essentially the same as that of the single- hole case. The dome structure is very 
pronounced. Further, our results shows that the amplitude of the background deviation, An = ((71.(0,0)) — { n (-n.-n)))i 
is roughly the same as (An^ hole + An^ hole ). This shows that the background behaviour is due to the single-hole 
excitations and is irrelevant to the physics of the bound state. We will provide a support to this in the next Sections. 

In the next section we will show that the important EMDF data are those along the AF Brillouin zone boundary. 
These data are practically unaffected by the kinematic form factor effect, so they can be used to draw conclusions 
on the internal structure of the bound state in k-space. One can interpret the EMDF at these points as the half- 
filled EMDF suppressed by the hole-occupation number. The hole weight at the single-hole ground state momentum 
(tt/2, 7r/2) is surprisingly small — (rik) deviates from the half-filled value of i by only 1%. This is the consequence of 
the d x 2_ y 2 symmetry which restrict the hole weight to be zero at these points. Another interesting feature is that the 
hole occupation at the (37r/4,7r/4) point is higher than that at (7r,0) (Fig. ||(a)). It is worth noting that at smaller 
J/t the hole occupation numbers at these points are very similar and their absolute values are larger. As it follows 
from the discussion in the next sections, these facts indicate the presence and importance of higher harmonics in the 
bound state at smaller J/t, because the "bare" first d-wave harmonic (cos(fca;) — cos(fc y )) will always give a larger hole 
weight at (tt, 0) than at (37r/4, 7r/4). 

The available clusters allow us to perform FSS for six of the nine inequivalent k-points of the 32-site cluster. Results 
for four of them at J/t = 0.3 and J/t = 1.0 are presented in Figs. ^|(a),(b). They all show the anticipated 1/N scaling. 



Note that a similar scaling plot at (ir/2, tt/2) is not successful because | (dn^) | is too small. Figure |10| shows the scaling 
of the EMDF at (7r,0). If we discard the 16-site data by arguing that they are spoiled by the artificial degeneracy, 
one can clearly see the 1/iV scaling at J/t = 1.0. In contrast to this, |(^(jr,o))l a ^ J/t = 0-3 does not show the same 
the 1/N scaling. We attribute these different behaviours to the different sizes of the bound states. The J/t =1.0 
bound state is small. Therefore, it has to scale as 1/N even when N is not too large. The J/t — 0.3 bound state is 
relatively large. An increase in the cluster size redistributes the hole weight among the new harmonics which become 
available on larger systems. The EMDF at those points not along the AFBZ boundary (Figs. ^|(a),(b)) mostly result 
from kinematic effects which are saturated at shorter distances. Therefore, they do not depend much on the details 
of the bound state structure. 
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FIG. 9. 1/N scaling of the |(<5rik)| in the two-hole ground state, for (a) k = (n,n) at J/t — 0.3 (filled circles) and J/t = 1.0 
(empty circles), and (0,0) at J/t = 0.3 (filled diamonds) and J/t = 1.0 (empty diamonds), and (b) k = (tx,-k/2) at J/t = 0.3 
(filled circles) and J/t = 1.0 (empty circles), and (7r/2, 0) at J/t — 0.3 (filled diamonds) and J/t — 1.0 (empty diamonds). 




FIG. 10. |(<5nk}| in the two-hole ground state vs. 1/N for k = (%, 0) at J/t — 0.3 (filled circles) and J/t — 1.0 (empty circles). 
Solid line shows 1/iV scaling for the J/t — 1.0 data if the 16-site cluster result is ignored. 
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IV. ANALYTICAL RESULTS 



Studies of the doped t — J model via the ED technique provide important information on effective quasiparticle 
theories. However, these same numerical results also posed some problems and made questionable the relation of these 
analytical studies to the problem of the "finite doping of the finite system" . For example, the EMDF for the ground 
states of the different number of holes and the pair correlation function for the two holes doped into the system were 
intensively studied numerically (see Sees. II. A, III). It turned out that the results for these quantities were found to 
be in the contradiction with some expectations. EMDF, which was naively expected to show something like "hole 
pockets" , or simply hole-rich and electron-rich regions in k-space, demonstrates a dramatic deviation of this quantity 
for the doped clusters from the half-filled (no holes) case with the strong variation across the whole Brillouin zone. 
Moreover, there is a strong doping dependence of these results. Data for the two-hole ground state differ significantly 
from the single- hole ones. More surprisingly, the overall shape of the EMDF reminds one of the free electrons with a 
nearest-neighbour hopping band. This was the reason for the conjecture that the t — J model already at rather low 
doping concentration undergoes a transition to the free-electron physics and shows a "large" Fermi surface |57| ], Also, 
the hole-hole correlation function for the (i-wave bound state show the largest weight of the holes in the configurations 
which should be forbidden by the d-wave symmetry (the so called "-s/2-paradox"). 

In this situation, a physical explanation of such puzzling behaviour of the correlation functions together with an 
analytical picture would be highly desirable. 

A qualitative understanding of these effects in the context of the spin-polaron physics has been achieved in the works 
by Eder and Becker |l^] , and Eder and Wrobel (l(| , wherein the authors showed that the t — J model quasiparticles 
will favour qualitatively the same EMDF as the ones found in the numerical calculations. Using rather general 
arguments, they demonstrated that the "large Fermi-surface" is a consequence of simple sum rules and a minimum of 
the total energy, and it is completely irrelevant to the problem of the real Fermi-surface identification. The main idea 
of these works is that the hole-pockets should be attributed to the quasiparticles, not to the "bare" holes. Since the 
renormalization is strong only a relatively small part of the polaron can be visualized in k-space as a fermion having 
the certain momentum. The "dressed" part of the spin-polaron is responsible for the background in the EMDF, which 
is spread over the entire BZ. More specifically, the EMDF does not only measure the lack of the electrons due to the 
centre of the polaron, but it also keeps track of the hole distribution inside the polaron. Similar physics has been 
discussed recently in Ref. . 

Within the same theoretical framework, i.e. a variational string approach, we mention that the pairing problem 
for two holes has been considered elsewhere JItJ and the source of the large probability of finding holes in the ground 
state along the diagonal of an elementary square can be explained by the large weight of the "holex (hole+1 spin 
flip)" combination in the two-hole <i-wave bound state wave function. Qualitative discussion of the same physics has 
been done recently in Ref. ]79[ |. 

In what follows we will show how the qualitative picture drawn in Refs. ||l5|-|l7f, which gives a basic understanding 
of the numerical data, can be reproduced using simple ansatzes for the spin polarons and their bound state. Then, 
the CT approach is used to derive analytical expressions which are able to explain quantitatively most of the 1- and 
2-holc ED data for the ground states described in Sec. II, and earlier in the literature. 



A. Qualitative analysis using a simplified model: 

We begin our analytical calculations of two holes described by the t — J model by considering a simplified treatment 
of two holes moving in an AFM Ising background. We then evaluate the EMDF and C (r) using this simplified model 
in order to see what kind of behaviour one can expect in the ground state of spin polarons. This work is instructive, 
and helps in understanding this problem. So we present these preliminary results first. 

Consider the EMDF. In a system without holes, one finds that (n^ a ) =1/2 everywhere in the full Brillouin 
zone. This is a consequence of the purely local character of the electronic states and Sf ot = 0. By definition 
(nko-) = jf e * kri3 (cLcj CT ) = Y,i(ctaCi<r) + W J2i,d^o e% *' d l (4a^i+du) • The second term is zero for the half-filled 
case and the first term yields (n^) — (rikj,) = 1/2. An informative result which follows from the second term 
is that hole doping makes the matrix elements between different "strings" of the polaron wave function nonzero, 
and accompanied by the phase factor e lk ' d , where |d| is the difference between the lengths of the strings. For 
example, the matrix element between the bare component Ciie lkri |0)) and the one spin-flip string component 

(Si is 5 , J + Ci+5te lkri |0)^ is proportional to X)a e * k ^ ~ 7k, which is asymmetric with respect to the transformation 
k — * k+(7r, 7r), 7k = — 7k+Q- Thus all odd-distance matrix elements are responsible for the antisymmetric contribution 
to (nk), and this asymmetry makes (rik) resemble the shape of a large Fermi surface. Note, that this unusual effect is 
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closely related to the localized character of the electronic states and the spin polaron nature of the carriers. Recently, 
a similar asymmetry observed in the ARPES data of an undoped (SrgCuOgClg) |^0| and doped JlJ] AFM insulators 
has been successfully explained using essentially the same ideas. 

Hole excitations near half-filling (when long-ranged AFM order is present) are most concisely explained using 
the spinless-hole, Schwinger-boson representation for the constrained fermion operators. Thus it is necessary to 
express (rik) and C'(r) in terms of averages of combinations of the hole and magnon operators. The essence of this 
representation is the following. The creation of a hole (annihilation of an electron) at site i in sublattice A = {f} 
(with the main direction of the spins being up) is achieved by operating on the ground state. Thus, c$f ~ h\. The 
action of the same operator on site j in sublattice B — {J,} is non-zero only if the spin is in the "wrong" direction 
(|). Therefore, creation of a hole is accompanied by the annihilation of a spin excitation: Cj-f = h^Sj ~ ^] a j- Thus, 

(4t^t) = (^A,ih A j(l - a' Aij aA,j))Si,A^j,A + {hA,ih Bj a B ,j)5i,A5j,B 

+ {hB,ta B , j h\ j )8 lt B5j,A + (hB,ih B] a Bj a B j)S lyB Sj, B (6) 

The above will be suffice for the description in this paper — for an advanced and detailed discussion of this represen- 
tation we refer the readers to Ref . (8lJ . 

First we examine a simple ansatz for the single-hole ground state wave function [p7|,[4l| 



i 1 ) = Vl^ |0) = V^ 



ah B V + 4/3 7P-q/4,P-q a B,q 



10) , (7) 



where a 2 + 4(3 2 = 1, and, as noted in Ref. the sign of the term linear in 7k is found from minimizing the kinetic 
energy. (Note that the origin of the hole is in sublattice B, so the total spin of the system is Sf ot = 1/2.) Hereafter, 
q € ABZ 

2 N/2 

Yl = Jq X] ' and { h A(B)k,h\ {B)k ,} = [a A ( B )k, a A (B)k>] = 5 W ■ N/2 . (8) 

q q n ^n— 1 

Minimal algebra for the EMDF and C(r) yields 

(n kT )A~ — - {h\ M h AM ) +XX( a kq as >q> ~ ( a iq a ^^q)) _ (( h A,k ^2 h B,k+i a B,q) + H - C -)' 



q 



(n k i)N~ y - (h BM h BM ) +^(( a iq a ^,q) _ ( a kq as <q)) ~ (( ft kk X] h A,k+n a A,q) +H.c), 

q q 

^"jSj^E^SVil.r). (9) 
id 

where = h\hi is the hole number operator. The physical meanings of the terms in (nk) are apparent. The number 
of electrons with spin up and momentum k is reduced by the amount of holes having the same momentum and by the 
spin flips in sublattice A. It is increased by the number of spin excitations in sublattice B. The last term is not zero 
between different components of the spin-polaron wave function, reflecting the inner structure of this quasiparticlc, 
or the kinematic "form- factor" . Alternatively, according to Ref. |19] , it reflects "the fast movement of the hole inside 
the bag" . 

Using |1) from Eq. (@), Eq. (|) give 

(n kT ) ~ \ + 1 (-16/3 2 7 2 + 4/3 2 + 8 | a ^| 7k ) f (10) 

These simple expressions already contain significant qualitative information about the EMDF for the single-hole 
ground state. There is a "dip" in (rikj.) at k = P with weight equals to one-half of the quasiparticle residue, 
corresponding to the centre of the polaron. There is also a constant positive (negative) shift in (n^) ((«k|)) due 
to spin excitations in sublattice B (see Eq. ([?])). Although (n^) does not have any "dips", it does have two other 
features. One is due to the hole distribution in the "dressed" part of the polaron (~ 7^), and the other is due to 
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the "interstring" matrix elements (~ 7k)- The absence of the "interstring" terms in (nkj,) in Eq. ( |10| ) is due to the 
approximation made in the above ansatz viz. Eq. (Q), namely the elimination of the longer strings which are necessary 
to produce the "dome" structure of (nkj.)- This explains the smaller amplitude and stronger t/J dependence of the 



difference An^ 



than those of An^ for the single-hole GS. 



Thus, most features of the single-hole EMDF data reported in Sec. Ill can be understood using this simplified 
model. According to Eq. (|Io| ) all (Sn^) scale as 1/N except for the dip which scales as C + a/N, a result which is 
employed in the FSS analysis of the ED numerical work. 

In order to carry out similar analysis for the two hole case, one has to solve Schrodinger equation for the bound 
state problem. Instead of doing this we simply propose the nearest-neighbour bound state wave function having 
d-wave symmetry and S% ot = based on the expectation that two static holes attract each other through the "sharing 
common link" effect (viz., H int = — J/2 n^rij): 



|2) 



v^£ A P^>Uio> 



p 

N 



p-q a B,q 



(11) 



16/3 2 7p-q7-p-q'/43,p-q / 4,-p-q a A,q a kq' 



|0> , 



with A£ 



cos(p x ) — cos(p y )) ensuring the centres of the polarons are at the nearest-neighbour sites. 
The EMDF calculation using |2) in Eq. (O) yields 



(r»k) = (n kT ) - (n kl ) ~ £ + 1 (-« 2 (A^) 2 - 16/? 2 7 £ + 8|a/3| 7k ) , 



(12) 



where the terms inside the bracket are simply the sum of the 1/N terms in the (nkf) and (riki) expressions of Eq. (|Io| ) 
for the single-hole case, and the "dip" structure is replaced by the probability of finding a "bare" hole with momentum 
k in the bound state. This explains the observation mentioned in Sec. Ill that the quantity An = nm.o) — n (Tr,Tr) 
for the two-hole case is roughly the same as An-[ + An^ for the single-hole case. It is interesting to note that since 
the k-dependence of the terms from the dressed part of the polaron and from the "interstring" processes vanish at 
the boundary of the magnetic Brillouin zone (7k = 0), features of the EMDF along this line are not disguised by 
kinematic effects. Thus, one can directly observe the structure of the bound state wave function A^ from the (rik) 
data at these points. In particular, the (tt/2, ir/2) point has to have zero hole weight due to the cZ-wave symmetry of 
the bound state. For the particular form of Ap we have chosen, the maximum of the hole weight (minimum in (rik)) 
will be at (jr, 0) point. 

The hole-hole correlation function on different clusters consistently shows a maximal probability for states in which 
the holes are along the diagonal of an elementary square, i.e., they prefer to be at a distance \l~2a from one another, 
where a is the lattice constant. (At first glance, such a configuration should be forbidden by the d x 2 _ y 2 -symmetry 
of the state. One way to resolve this paradox, as suggested by Poilblanc @|, is to introduce modified creation pair 
operators '^jixiy'-'ii-cfy) ^° ^ ne "b are " ^l^lixfy) P a ^ r °P era t° r - It is clear that the bound state wave function Eq. 
( pi] ) includes such combinations naturally.) Calculation of C(r) of Eq. (^) in the ground state given by Eq. ( |ll| ) gives 



C(l) = a 4 /4 + 9/3 4 /4, C(V2) = a 2 (3 2 , C(2) = a 2 (3 2 /2, C(Vb) = 3/3 4 /4, C(3) = /3 4 /4. 



(13) 



3 the weights accumulated in the 
y 2 ~ A(3 2 MM. This gives C(l) 2 



'bare" and "1-string" parts of the polaron 
C(y/2), in qualitative agreement with the 



For the physical range of t/J ~ 2 
wave function are almost identical, 
numerical results. 

Thus, one can conclude that our simple considerations of one and two holes in a system of Ising spins, based on 
a simplified spin-polaron picture, already shows qualitative agreement with the numerical data. The treatment of 
the realistic system with a Neel spin background requires a proper account of the spin fluctuations, the long-range 
dynamics of the system, and multiple spin excitations (longer "strings"). 



B. CT approach 
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1. CT Hamiltonian: 



The t — J model Hamiltonian ([!]) can be rewritten using the spinless-fermion representation for the constrained 
fermion operators and Holstein-Primakoff (l^jlj] or Dyson-Maleev representation for the spin operators. These 
formalisms have been shown to be adequate in treating the nonlinear feature of the kinetic energy term of Eq. (lj) 
properly. Subsequent diagonalization of the spin part of the Hamiltonian by the Bogoliubov transformation naturally 
includes spin-fluctuations in the ground state. 

The essential part of the t — J Hamiltonian rewritten in this way looks like the electron-phonon Hamiltonian for 
the "usual" polaron problem with an additional direct fermion-fermion interaction term: 

Ht-j =s w q a i Q q + X! { M ^ h l-^ a l + H - c -) + AH > ( 14 ) 

AH = -2.7(1 - 2A) Y, 7q4- q 4' +q ^k , 

k,k',q 

where h\h), c^(a), are the spinless hole and magnon operators, respectively, w q = 2J(1 — 7q) 1 ^ 2 is the spin-wave 
energy, Mk, q = 4i(7k- q E/ q + 7k Vq), C/ q , Vq are the Bogoliubov transformation parameters, 7k = (cos k x + cos k v )/2, 
A = J2q(Vq ~ 7q^q^q) = —0.08. Ai7 is an effective hole-hole attraction due to minimization of the number of broken 
AF bonds. Two important differences make the t — J version of the polaron problem much more difficult to study: (i) 
the absence of "bare" dispersion term of the hole p^] , and (ii) the essentially non-local character of the hole- magnon 
interaction, because each process of emitting (absorbing) a magnon is associated with an intersite hole hopping. 

The CT approach has been applied to the spin-polaron Hamiltonian of Eq. (|l4|) in Ref. jl9| . The generator of the 
CT was proposed to be in the form 

s = Y /a,, (4-,fck°4 - Hx -) > ( 15 ) 

k.q 

where /k is the parameter of the transformation, and in Ref. Jl9| /k was chosen to minimize the single hole energy, 
viz. 



S 



(i>)=° 



(16) 



The negligible role of the higher order hole-magnon vertices in the transformed Hamiltonian was demonstrated and 
it was argued that the initially strong hole-magnon interaction in Eq. (JL4|) is transfered mainly into a hole " dressing" 
and into the hole-hole interaction. Thus, for a wide region of tj J one can restrict one's considerations to the effective 
Hamiltonian 

H -ff = X + X w q« q «q + X Vw^-Ji'+Jivh* + ^ q M k ,q (ht_Ji±al + H.c.) (17) 

k q k,k',q k,q 

where and w q are the polaron and magnon energies respectively, Vk,k',q is the direct polaron-polaron interaction, 
A^k,q is the bare hole-magnon vertex, and Fk, q is the renormalization form-factor which is close to zero at large 
q, and is constant (~ 0.2 — 0.4) at small q. The last term, which corresponds to the interaction of the hole with 
long-range spin waves, has been left in the effective Hamiltonian in this form to account for the retardation effect 
in the polaron-polaron spin-wave exchange. Also, short-range spin-wave exchange has been converted to the direct 
polaron-polaron interaction. The polaron energy and the weights of the components of the polaron's wave function 
have been compared with the results of the other works, especially SCBA results, and very good agreements were 
found. Since the derivation of the single polaron energy and the polaron-polaron interaction in the framework of CT 
approach are the same, one can hope that the effective Hamiltonian of Eq. ([l7]) properly describes the interaction 
between the low-energy excitations of the t-J model. 

2. Our calculations using the CT approach 

We are interested in the ground state with total spin S^ ot = 1/2 (S^ ot = 0) for the single-hole (two-hole) case in 
an AF ordered system. Thus, it is necessary to use a two-sublattice representation for the fermions and bosons p9[. 
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In the two-sublattice representation there are two types of holes and magnons, both defined inside the first magnetic 
Brillouin zone, whereas in the one-sublattice representations holes and magnons are defined inside the full Brillouin 
zone. In the previous subsection we used the latter for the sake of simplifying notations. There is a simple relation 
between these two representations: 

h k = (/k + 9k)/ V2 , /ik+^vr) = (/k - <7k)/\/2 , (18) 

Oq = (a q + 0q)/ V2 , flq+^.Tr) = (Oq - /3q)/ ^2 , 

where / k and g k correspond to the hole excitations in the A and B sublattices respectively. a q and /3 q are the two 
branches of the Bogoliubov spin- wave excitations. 

The correlation functions (n k(T ) and C(r) expressed in terms of the averages of the hole and magnon operators are 

Kt)^ * f - (1 - SX){fth) E^+qSk+q)^ + J2 (^q/U ~ «Oq>) 

q q 

- ((/^.gk + q(/9 q C/q + a-q^q)> +H.C.), 

q 

(n n )N ~ y - (1- «5A)( 9 t 5k ) - ^</ 1 l +q /k + q)V q 2 + E «« q a q ) - 

q q 

- ((sEW^+^q))+H.C.), 

q 

C(r) = ^ (r) X>N£*|i-i|,r>, ft = /(<?), when * e {A}({5}), (19) 

where = 5Z q = 0.19. The negligible contribution of the higher order terms (in the number of magnons) to 
(nk,<r) has been checked and these terms are omitted. 

It is interesting to compare these expressions with those for the Ising limit of the model given in Eq. (|l(]). The 
number of holes in sublattice A reducing the number of electrons with spin up is decreased by the spin fluctuations 
(<5A, first term), but due to the same effect the reduction in (n k -|-) can be done by the holes in the sublattice B (second 
term). The third terms take into account an imbalance of the number of spin excitations of different types. The last 
term is nonzero for the "interstring" processes (now "strings" are just the components of the wave function with the 
spin excitations). 

Using the CT generator of Eq. (p"5|) one obtains the wave function of the spin polaron {Sf ot = +1/2) 



|1> = V^5pI°> = \[^ eS 9>) = [ap ffP + £ 6p,qJp- q /?q + E CP,q, q ' <7 P _ q _ q ,/?H + 



10} ■ (20) 



Here, ap = Zp < 1 is the quasiparticle residue. An explicit expression for the exact spin-polaron wave function within 
the SCBA was written in Ref. p4j in the same form. The ground state momenta for the spin-polaron in the pure 
t — J model are ±(±7r/2, ir/2). 

Then, using Eq. (p0|) the single-hole EMDF is found to be given 



(n kT ) ~ | + — ( -(1 - SX)b p> p_ k - fl^P-k - E C P,q',P-k-q-q'Kj 2 

q,q' 



2 N 

+ E fo P,q _ 2(a P &p,p_ k J7p_k + frp,p- k E c P . P _ k!q V q ) J , (21) 
q q ' 

M ^k,p(l - 5X)al + 1 (-(1 - SX) 4,q,P-k-q - E ^.P-k-q^ 

q q 

— E ^P,q _ ^ E kp,P-k-q C P,P-k-q,q^q ) • 

q q 

As we will show below, these expressions give good quantitative agreements with numerical data. As before, (n k j.) 
shows a "dip" at k = P with a weight proportional to of the quasiparticle residue due to the centre of the polaron. 
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A constant positive (negative) shift due to different amount of spin excitations (fourth term) is also present in (nkf) 
((riki))- The first three terms in (n^) and the second and third terms in (rikj.) reflect the hole distribution inside the 
polaron, whereas the last two terms in (n^) and the last term in (rikj.) are from the "interstring" matrix elements. As 
before, they are odd with respect to the transformation k — * k + Q and lead to the formation of the "dome" structure 
in the EMDF. As we noted, the asymmetric term in (rikj,) comes from the matrix element between the second and 
third components of the wave function of Eq. (^o|) . We restrict ourselves to the first three components of Eq. (|2(i| ) 
because for J/t = 0.3 they give about 98% of the norm of the wave function. (We note that in the SCBA approach 
the same approximation gives about 92% of the norm ]32]| ). 

Formally, Eqs. @ give a 1/N scaling for (Srika) a t every k-point except for the "dip" in (<5nkj,) at k = P . It 
fails at the point k = P — Qaf where some of the terms in Eq. ( ^l| ) are singular. The reason for these singularities 
is a peculiarity of the spin-polaron ground state and the AF long-range order. The dressing of the hole in the Neel 
background involves an infinite amount of zero-energy q = Q spin excitations (whose total contribution to the hole 
weight is finite and small due to the diminishing magnon density of states). Since the EMDF probes the inner structure 
of the ground state it is actually measuring this singular probability of the virtual emission of a zero-energy magnon 
(Q) by the hole (P) if k is equal to P — Q J8j|. This leads to singularities of different types for (rikt) and (n^i): 

Kr> - L(k-(P- Q )r M ~ a7 ln(w(k - (P - Q)) < (22) 

where u(k) is the magnon energy. For the finite system the magnon spectrum has the finite energy gap at Qaf which 
scales as [ p4| , [S5[ 

AE = J*'(l-±^' + ..\, (23) 



p s N V Ps 4tt 

where c ~ 1.67 and p s ~ 0.175 p6[ are the spin-wave velocity and spin stiffness, respectively. This result gives 
"antidips" reported in Sec. Ill with the following scaling laws: 



-t(p-Q)>^ + §-£o T -^L + ^, <-i(p-Q)>4-^- 



2 N Jn N 2 N N 



where all constants are positive. Cj, C| are from the "regular" part of Eq. ( |2l[ ). An interesting result shown in 



Eq. ( |24| ) is that the "antidip" in (rtj(k)) is predicted to survive in the thermodynamic limit: 

(n T (-7r/2, -vr/2)) ~ Z P p s /c 2 ~ 0.07 • Z P , (25) 

whereas all other features except for the dip at P in (nj(k)) ~ Zp will disappear. One can see from Eq. (|2~i| ) that 
the scaling laws of the antidips are quite complicated. For a system as small as N = 32, terms of different order in 
N have similar amplitudes. For example, By\/y&, ~ 0.5-Bot- This makes the FSS for the "antidips" complicated, 
especially when only two of available clusters (16 and 32) possess this k-point. An additional complication comes 
from the fact that the gap AE Eq. ( p3| ) is calculated for the system without holes and the influence of the latter on it 
is not known. In the subsequent calculation of the EMDF for "antidips" we modify the magnon spectrum employed 
in Eq. (|l]) in a way that it has a gap AE at q = Q point |87j . 

The two- hole problem has been considered using the Hamiltonian of Eq. ( |l7|) in Ref. fl8[| . A bound state with 
d x 2_ y 2 symmetry was found for < t/J < 5. The wave function of the d-wave bound state with total momentum 
P = can be written in the terms of creation operators of polarons of Eq. (|2p|), 



12) = |^ =0 ) = ^1 E A p/p.9-pl°> ( 26 ) 

P 

OO OO 

A p = E C 2n -i,2m{cos([2n-l}p x + 2mpy) - cos(2mp x + [2n-l]p v )} 

n—1 m=— oo 

= Ci,o {cos(p z ) - cos(p y )} + C 3! o {cos(3p x ) - cos(3p y )} 
+ {cos(p x ± 2p y ) - coa(p y ± 2p x )} + . . . , C 2n -i,2m = C 2n -i-2 m 

where Ap is the solution of an analog of the Schrodinger equation for the two-body problem. The form of Ap ensures 
the d-wave symmetry of the state and that the centres of the polarons are always on different sublattices, which in 
turn guarantees S z = 0. Ap has a more general form than the simple "nearest-neighbour c?-wave" in the simplified 
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example of Eq. (jllj). It has higher harmonics, which have substantial weight for realistic t/J. In what follows we 
show that our comparison leads us to the conclusion that the large higher harmonics of Ap play an important role in 
determining the behaviour of C(r) and (rik)- 

Note that for the representative value J/t = 0.3 about 42% of the polarons in the bound state are located at the 
nearest-neighbour sites and less than 2% are farther than 7 lattice spaces. An interesting feature of this distribution 
is that the probability in finding the second polaron at a certain distance from the first falls off slower along the x 
and y directions. Thus, the weight of the (3,0) (cos(3p x ) — cos(3pj,)) harmonic of Ap is rather large (20%), whereas 
the weight of (1,2) (cos^ ± 2p y ) — cos(2p x ±p y )) component is less than 5%. 

Finally, relating the previous forms of the wave function (in terms of creation operators of holes and magnons) Eq. 
(p6|) becomes 



|2> 



At(p)+At(p)+At(p) 



|0> 



(27) 



p+q 



ft 



(2) 



,t „t 



P:q,q '3p+q'5-p-q-q' 



A|(P) = E 

q,q' 



n W f t „t ,r|W f t„t , C W J ft 

w p,q,q'Jp — q— q'^— p ~ — p,qf p&— p— q— q' ~ p,q,q' ^p~q^ — p — q' 



where the subscripts n of A's indicate the number of magnons in the corresponding component of the wave function. 

Results of the EMDF and C(r) in Eq. ([l9|) for the ground state of Eq. ( p7| ) are given in full detail in the Appendix. 
In the next section we use these expressions to compare this theory with the numerical data discussed earlier in this 
paper. 



V. COMPARISON OF NUMERICAL AND ANALYTICAL RESULTS 



This section summarizes the comparison of our numerical ED data with the analytical results obtained from the CT 
approach. We focus on the EMDF for one and two holes, the binding energy, and the hole-hole correlation function 
for two holes. These provide a representative juxtaposition of results obtained from these two techniques, and probe 
in detail the correlations found in the ground states. 

Figures mjl2|Jl^ , [l4l show our analytical results for the single-hole EMDF (Eq. (^l|)) together with the 32-site ED 
data. Solid lines are guides to the eye. The agreement is excellent for both spin directions. The differences between 
the (Jrikt) numerical and analytical data at (0,0) and (tt,tv) can be attributed to the fact that the CT quasiparticle 
residue = a£ at these points is larger than the "exact" values (e.g., SCBA). As one can see from Eq. (21), this 
leads to lower values of (Sn^). The agreement of the (Sn^i) quantities away from k = P is better because the role 
of the background does not depend on a^. The "antidips" in the analytical results are marked by the cross notifying 
that these points were calculated from Eq. (^) using finite gap value in the magnon spectrum \ pT\ . 

As we discussed in Sec. IV, the internal structure of the spin polaron is made evident in the EMDF through the 
"normal" and "interstring" terms. The normal terms reflect the distribution of the hole inside of the spin-polaron 
wave function, viz. strings of different length. Interstring matrix elements arc non-zero for (#«k) due to the specific 
structure of the spin polaron. They make a contribution to the EMDF Eqs.([l0|),(^l|) which is asymmetric under the 
transformation k — * k + Q. These asymmetric terms are responsible for the dome shape of (Sn^), as was proposed 
in Ref. thus showing that it is not related to a Fermi surface signature. 
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(-71,-71) (0,0) 



(7t,7t) (7C,7C) 



(— 7t,7C) (— 7t,7t) (0,0) (7E,— 7U) 



FIG. 11. Comparison of the numerical (open circles) and analytical (filled diamonds) results for (Sn^i) in the single-hole 
ground state at J/t = 0.3 along the lines shown in the inset ((— n, —n) — > (n, n) — > (— n, n) — * (n, — tt)). Solid lines are a guides 
to the eye. 




(7C.0) (71/2,71/2) (0.7C) (0,7C) (-71/2,71/2) (-7C.0) (-7i;0) (-71/2,-71/2) (0,-71) 



FIG. 12. The same as in Fig. [H] along the lines ((tt, 0) -> (0, tt) -> (-tt, 0) -f (0, -tt)). 



19 



0.02 



eg 
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(0,0) 
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(c) 



(7C,7C) (7C7C) 



(-7t,7C) (-7C.7C) 



(0,0) 



(7t;-7T) 



FIG. 13. Comparison of the numerical (open circles) and analytical (filled diamonds) results for (Sn^i) in the single-hole 
ground state at J/t = 0.3 along the lines shown in the inset ((— n, —n) — > (n, n) — > (— n, n) — > (n, — 7r)). Solid lines are a guides 
to the eye. 




FIG. 14. The same as in Fig. |l| along the lines ((tt, 0) -> (0, tt) -> (— k, 0) -> (0, -tt)). 

The J/t dependence of the single-hole EMDF data has been extensively studied in Ref. |38| for smaller systems. 
As we already noted in Sec. IV, Eq. (|2l|) naturally describes the results of these studies and of the observations made 
in Sec. III. A. To be specific, the depth of the "dip" is proportional to Zp, and thus (Snp^) must follow the J/t 
dependence of Zp (result of Ref. |[|]). Also, the background must be getting weaker for larger J/t because less hole 
weight is accumulated in the string cloud. The smaller values and stronger J/t dependence of A?t,| and An ant i,i arc 
due to the second and third components of the spin-polaron wave function involved in the formation of the (Sn^^) 
background, which are more sensitive to J/t. 

Now we consider two hole results, beginning with the binding energy. This quantity for cZ-wave bound states was 
obtained numerically and analytically at a variety of J/t. For J/t = 0.3 it was found that Ef D = -0.05t and 
E^ T = — 0.02t. As discussed before, the absence of a simple scaling law for Ej, does not allow one to produce a 
reliable estimate of its thermodynamic limit at small J/t. Nevertheless, we believe that the close agreement of the 
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energies supports the idea that the systems under study represent the same physics. A simple 1/N FSS for the larger 
J/t = 1.0, where the size of the bound state is smaller and it is hoped that such a FSS will be more credible, gives a 
thermodynamic value of E^ D ~ — 0.32i (FigJ|(b)), which is very close to the theoretical result Eff t — — 0.38£. 

Figure 0(a) shows our results for the hole-hole correlation function C(r) for two holes in the d-wave bound state, 
for J/t = 0.3. This expectation value is calculated firstly for a bulk lattice, and then mapped onto the equivalent sites 
of a 32-site cluster with periodic boundary conditions. This enforces that the analytical work approximates some of 
the finite-size effects of our ED numerics, and facilitates a more natural comparison between the two. Very similar 
trends are found in both results, with the correlation function decreasing quite similarly with the distance. 

Both numerical and analytical results data show that about 45% of the time the holes prefer to stay at the nearest 
and next-nearest neighbour distances. That our analytical work produces such behaviour is not inconsistent with our 
statement regarding the form of Ap Eq. (|2^) : the centres of the polarons are indeed restricted to be on opposite sub- 
lattices, but the holes are almost equally distributed on both sublattices, with the maximum probability of separation 
being at y/2. In fact, our analysis of the harmonics in Ap given in Table I shows that about 40% of the polarons in 

the bound state are separated by one lattice constant. Thus, the peak of C(r) at r — v2 arises from the components 
of the wave function with "strings" of length one. Clearly, the spin-polaron picture provides a natural explanation 
for the -\/2-paradox found here and in earlier numerical studies. 

Table I. Amplitudes (C2 n -i,2m) and weights (Cf„_i 2m) of the harmonics in the spin-polaron bound state wave function 
Eq. (p6[). Weights are directly related to the polaron-polaron spatial distribution function: P(r u > ) — -^-jC^, , Zu> — 4; 8 is the 
coordination number. / and/' belong to different sublattices 



(2n- 1,2m) 


G^n— 1,2m 


C2n-1.2m, % 


(1,0) 


0.642 


41.3 


(1,2) + (1,-2) 


0.215 


4.6 


(3,0) 


0.444 


19.7 


(3, 2) + (3, -2) 


0.320 


10.2 


(5,0) 


0.240 


5.8 



We believe that there are two reasons for the analytical (7(1) being slightly larger than C(V2). A treatment of the 
t- J model based on the spinless hole representation involves some unphysical states with the hole and spin excitations 
being at the same site. The number of processes leading to such states increases when the polarons are close to each 
other and hence (7(1) grows. Secondly, the CT approach slightly overestimates bare hole weight. 

An additional maximum in our analytical (7(r) at r = 3 is closely connected to the second important harmonic in 
the CT d x 2_ y 2 bound state ~ (cos(3pa;) — cos^p^)). We cannot explain the absence of such a peak in the ED results 
— for example, we have been unable to estimate the finite size effect on individual harmonics in the two-hole wave 
function. 

Figure ^(b) demonstrates a better agreement at J/t = 0.8 when the size of the bound state is small. In this case 
the correlation fall rapidly with distance and thus the "bulk to cluster" mapping does not alter the analytical data. 
Therefore, in the large J/t limit we find the expected result that a spin-polaron approach adequately describes the 
physics. 

The EMDF for two holes shows an equally satisfactory comparison, as seen in Figure |l5|. The behaviour of (Snu) 
involves the combined effects of (i) the internal structure of the polarons, (ii) the <i-wave symmetry bound state, and 
(iii) higher harmonics in the bound-state wave function. We now elaborate on these features. 

As we argued in Sec. III.B, the quantity An 2hole = «Ti( ,o)) ~ (^(tt.tt))) - (An 1 1 hole + &n\ hole ) shows that the overall 
background deviation is mainly irrelevant to the bound-state k-structure. The worst agreement between ED and CT 
analytics for the (0, 0) and (-7T, it) points (Fig. |l5|) is again due to the Zk =0 problem within the CT approach. 

Next we focus on the features along the ABZ boundary which, as we discussed before, are not disguised by kinematic 
effects (all asymmetric and most of the "normal" terms are zero on this line) and can be directly related to the form 
of Ap in Eq. (p6|). Our analytical and numerical (Srik) results have a local maximum at k = (tt/2, 7r/2), and have 
minima at (37r/4, 7r/4) and (7r/4, 37t/4). The first feature can be explained by the d-wave symmetry of the bound 
state. The EMDF is reduced from its half-filled value of 1/2 when holes occupy that momentum state. However, as 
shown in Eq. (|l2]), the EMDF consists of terms proportional to (A^) 2 , which is identically zero at (tt/2, 7r/2). Thus, 
(5nk) must show a local maximum (=0) at this wave vector, so one cannot find any direct remnant of hole pockets 
for d x 2_ y 2-w&ve symmetry bound state. 

A second feature that we observe in both analytical and numerical results, viz. the minimum along the ABZ 
boundary between (w, 0) and (tt/2, tt/2), can be related to the particular form of Ap. Analytically, this quantity has 
large and apparently important higher harmonics (see Table I and Eq. (p6h), and it is the competition between the 



21 



different harmonics that produces the maximum hole number between (7r/2, tt/2) and (37r/4, 7r/4). Our analytical 
work shows that the hole number is actually maximized along the ABZ boundary very close to (7r/2, 7r/2), roughly at 
(0.457T, 0.557r). It is unclear if experiments could resolve this feature. 



0.04 




-0.12 1 1 ' ' < < < < < 1 ' — 1 

(0,0) (7E'2,7Ii'2) (Tt,Jt) (71,71) (0.7C) (0,7t) (7172,71/2) (7C.0) (7C.0) (0,0) 



FIG. 15. Comparison of the numerical (open circles) and analytical (filled diamonds) results for (Sn^i) in the two-hole ground 
state at J/t — 0.3 along the lines shown in the inset ((0, 0) — ► (tv, n) — > (0, tt) — ► (tt, 0) — > (0, 0)). Solid lines are a guides to the 
eye. 

VI. CONCLUSIONS 

Summarizing, we have presented new ED numerical data for up to two holes in the t-J model for the largest 
cluster for which such calculations can be completed presently. Then, we compared these results with new analytical 
expressions based on the canonical transformation approach to the t-J model. We find good agreement for the binding 
energy, the EMDF for one and two holes, and the hole-hole spatial correlation function. We consider this to lend 
strong support to the validity of the quasiparticle Hamiltonian derived in Ref . [ fl9| , thus supporting the contention 
that the spin-polaron description of the quasiparticles in the t-J model is correct at least at low hole concentration. 

Certain characteristics in the correlation functions we studied are direct consequences of some features of the 
corresponding ground state wave functions. For example, the dip in the single-hole spin J, EMDF is related to the 
centre of the spin polaron, whereas the dome structure of the 1- and 2-hole EMDF is due to the inter-string matrix 
elements of (nk CT )- The large correlation between the holes in the c£-wave bound state at a distance of \[2 is due to 
the significant weight of the shortest string in the spin-polaron wave function. Analytical results for the FSS of the 
EMDF show a 1/N scaling at almost all k-points except at the single-hole ground state P (Eqs. ([Tol) , (|l^) , ) , in 
agreement with what is shown in Sec. Ill for the ED data. Those k-points which are influenced by the long-range 
physics of the system are shown to have more complicated scaling laws — see Eqs. (^2|) , (|24|) . The role of the higher 
harmonics and the effect of the size of the bound state on the EMDF and the hole-hole correlation function for the 
two-hole problem are also discussed. 
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APPENDIX A: EMDF AND C(R) FOR TWO-HOLE CASE 



Amplitudes of the components of the two-hole wave function (§7|) A^ n \ flW, can be expressed through the 
a,b,c components of the single-hole wave function fed): 



) — fl k J ^k,q — 26k,q6-k+q,q , ^q.q' — c k,q,q'C-k+q+q',q',q 

= -a k 6_ k ,q , ^q.q' = ^k,q'fok-q',q6-k+q',q' (Al) 

G k,q,q' = 2 akCk .q^' ' °k,q,q' = Ok,q0-k,q'- 

flk,^k,q, and Ck, q ,q' within the CT approach are the products of Mk, q 's (Q), /k's, and different integrals of their 
combinations as described in Ref. . /k is the transformation parameter obtained from Eq. ( |l6| ) . 

Using Eq. (|H| ) and the bound state wave function (^) one can get EMDF for the two-hole ground state 



/ N 1 1 

^^2 + N 



-(8n^ en ) + 2{5n% Ld ) 



(Sum =FAl + Y, FBl^ + ]T FC£ q , q , (A2) 
q q,q' 

(5< d ) = FB^U^FA^ - FA k £ V^FB k , 



,q 

q q 



^ FCk^q'Kl'-F-E'k+q+q'^q + ^ FBk,q' VqFCk+q' , -q, -q' 

q,q' q' q 



where 



= + £ A£ +q 4 2 j q , q + £ A^ +q+q ,4 3 j q+q , q , q , 



q.q' 



FB k: q - A^B^ q - A^ +q S^ q _ q + ^ (Ak + q+q'-Sk+q+q',q,q' ~ &t+qf B k+q> ) ( A3 ) 

q' 

JU M,q' — ^k+q+q'°k+q+q',q,q' + ^k°-k,q' : q ^k+q' °k+q' , -q ,-q' ■ 



C 00 (m ) + C 22 (ry ) = (2 1 n[n] | 2) , when i + j = 2n - 1 
C n (?Xj ) = ( 2 kfn||2>, when i+j = 2n 



Hole-hole correlation function (|19|) in the bound state (|27|) is given by 
C(r l3 ) 

C °"(^) = (E A k{4 1) coslkr^] + £ 4 2 q cos[(k - q)r 4i ] + £ cos[(k q q')r y ]l) ' 

^ k q q.q' ^ ' 

C l \ rij ) = £ A£A* WfiW flW (cos[(k - p)r«] - cos[(k + p + q)r«]) (A4) 
k,p q ^ 

+ 2S« ^ fig iq/ (cos[(k - p + q + q')r i3 -] - cos[(k + p q> y ]) | 
q' 

k,p q,q' ^ ^ 



2C £' >q C -iq',q cos[(k p + q + q')r<,-] - 2C< 2 q , q C« _ q , cos[(k + p q')r tf ] 



241q'^ q , q cos[(k + p + q')r y -]| 
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